We wish to create conservation-sensitive, climate-smart fisheries closures for the ABNJ of the Pacific Ocean using prioritizr.
Prioritizr (Hanson et al., 2018) is a spatial prioritization software (R package) that uses integer linear programming (ILP) allowing it to create exact solutions to large spatial planning problems faster than its heuristic counterparts (e.g. Marxan).
This file runs multiple scripts found in the main Github repository: tinbuenafe/PacificProj, where the codes are explained in greater detail.
Before running prioritizr, we first had to assemble the following data layers:
The study area has the following characteristics:
We ran the following code using polygons from rnaturalearth package and the EEZ .shp files from Marine Regions (Flanders Marine Institute, 2019) to create a Pacific-centered map of the global ABNJs.
The relevant scripts for creating the study area are labeled with a prefix of 01 in ~/scripts/.
# to create the Pacific-centered maps:
source("scripts/01a_StudyArea_CreatingPUs.R")
# this script uses the following functions (embedded actual script):
# function to create polygon of the ABNJ by masking the EEZs and the landmasses
source("scripts/study_area/fCreateMaskedPolygon.R")
# function to reproject these polygon to Robinson projection
source("scripts/study_area/fConvert2PacificRobinson.R")
# function to create the hexagonal planning units
source("scripts/study_area/fCreate_PlanningUnits.R")
This figure shows the Pacific-centered global map (using Robinson projection).
world_robinson <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificCenterLand.rds")
global_plot <- ggplot() +
geom_sf(data = world_robinson, colour = "grey20", fill="grey20", size = 0.1, show.legend = "line") +
theme_bw()
global_plot
We then created equal-area hexagonal planning units, or PUs, (0.5\(^\circ\) x 0.5\(^\circ\) at the equator) of the study area. This resulted in 31,917 PUs.
# define objects for plotting; objects are the result of running 01_StudyArea
pacific_robinson <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificCenterABNJ.rds")
PUsPac <- readRDS("outputs/01_StudyArea/01a_StudyArea/PacificABNJGrid_05deg.rds")
Bndry <- fCreateRobinsonBoundary(west = 78, east = 140, north = 51, south = 60) # coordinates in degrees
study_area <- ggplot() +
geom_sf(data = pacific_robinson, colour = NA, fill = NA, size = 0.2, show.legend = "line") +
geom_sf(data = world_robinson, color = "grey20", fill="grey20", size = 0.1, show.legend = "line") +
geom_sf(data = PUsPac, colour = "grey64", aes(fill = "ABNJ"), size = 0.1, show.legend = TRUE) +
scale_fill_manual(name = "Study Area", values = c("ABNJ" = "coral3")) +
coord_sf(xlim = c(st_bbox(Bndry)$xmin, st_bbox(Bndry)$xmax), # Set limits based on Bndry bbox
ylim = c(st_bbox(Bndry)$ymin, st_bbox(Bndry)$ymax),
expand = TRUE) +
theme_bw()
study_area
After demarcating the study area, we categorized the PUs into their Longhurst Provinces, .shp files from Marine Regions (Flanders Marine Institute, 2009), using the function fCategProv() found in this source code.
# function to assign the provinces defined in:
source("scripts/01b_StudyArea_fCategProv.R")
# running this function:
source("scripts/01c_StudyArea_fCategProvRun.R")
longhurst
The following were used in this study:
First, we compiled a list of bycatch reported in the Pacific ABNJ. Then, we extracted their AquaMaps distribution. AquaMaps (Kaschner et al., 2019) is a niche-based model that predicts the relative occurrences of species from their environmental envelopes or their environmental and oceanographic preferences (Kesner-Reyes et al., 2016).
We limited the bycatch included in this study to sea turtles, but this workflow can be easily replicated and expanded to include more species.
The following were the sea turtles reported as bycatch and intersected with the study area:
The relevant scripts for creating the bycatch distribution layer are labeled with a prefix of 02_RawAQM and 03_AQM in ~/scripts/.
The plots below show the distribution of these species using a probability threshold of 0.5 and the functions fAquaStart() and fBycatchDistMap().
# to extract the distribution of all species that are included in the AquaMaps dataset:
# the function is found in:
source("scripts/02a_RawAQM_fAquaStart.R")
# the function is ran in:
source("scripts/02b_RawAWM_fAquaStartRun.R")
# to create the distribution maps of the sea turtle species included in this study:
# the function is found in:
source("scripts/02c_RawAQM_fBycatchDistMap.R")
# the function is ran in:
source("scripts/02d_RawAQM_fBycatchDistMapRun.R")
Just like the bycatch layer, we compiled a list of commercially targeted large pelagic species managed in the Pacific from different sources. From this list, we have 4 species with available data on their global larval spawning grounds:
The models for these species are from a report by Mercer (2019) where he created Generalized Additive Models (GAMs) to predict spawning grounds. He used average seasonal abundance from a historical dataset of longline fisheries surveys (dated back to the 1950s), and oceanographic and environmental conditions from global databases.
We reran his models (see scripts labeled with a prefix 05a in ~/scripts/ for complete reruns) and used the medians of the predictions for each species as probability thresholds to create the distribution maps and the function fCommercialFeat().
# to intersect the results of the GAMs with the PUs
# the function is in:
source("scripts/05b_Commercial_fCommercialFeat.R")
# the function is ran at:
source("scripts/05c_Commercial_fCommercialFeatRun.R")
global_plots
Each province is expected to be biogeographically and ecologically unique. Including this layer in the spatial prioritization spreads the risk of exposure of populations to calamities, diseases or pollution. We joined the features’ layers with the provinces’ layer, creating multifaceted features, with each feature subtyped by provinces. This increased the number of features from 9 (5 from bycatch and 4 from large pelagics) to 66 features.
We used the function fProvIntersect().
# to intersect the conservation features with the provinces
source("scripts/03_Features_fProvIntersect.R")
# function ran with the bycatch layer:
source("scripts/03b_AQM_fProvIntersectRun.R")
# function ran with the commercially targeted species layer:
source("scripts/05e_Commercial_fProvIntersect.R")
This is an example of how the data looks like (so far) for bycatch:
head(bycatch)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -3779074 ymin: 168708.9 xmax: -3695824 ymax: 473118.3
## CRS: +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
## # A tibble: 6 x 5
## # Groups: province [1]
## province cellsID species geometry new_feature
## <chr> <int> <chr> <POLYGON [m]> <chr>
## 1 WARM 1 ITS-Mam-… ((-3751324 216773.5, -3779074 232795.… WARM_ITS-Ma…
## 2 WARM 2 ITS-Mam-… ((-3751324 312902.8, -3779074 328924.… WARM_ITS-Ma…
## 3 WARM 3 ITS-Mam-… ((-3751324 409032.1, -3779074 425053.… WARM_ITS-Ma…
## 4 WARM 4 ITS-Mam-… ((-3723574 168708.9, -3751324 184730.… WARM_ITS-Ma…
## 5 WARM 5 ITS-Mam-… ((-3723574 264838.2, -3751324 280859.… WARM_ITS-Ma…
## 6 WARM 6 ITS-Mam-… ((-3723574 360967.5, -3751324 376989,… WARM_ITS-Ma…
And for the large pelagics:
head(pelagic)
## Simple feature collection with 6 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -3779074 ymin: 216773.5 xmax: -3695824 ymax: 1482476
## CRS: +proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
## # A tibble: 6 x 5
## # Groups: province [2]
## province cellsID species geometry new_feature
## <chr> <int> <chr> <POLYGON [m]> <chr>
## 1 WARM 1 ALB ((-3751324 216773.5, -3779074 232795.1, … WARM_ALB
## 2 WARM 2 ALB ((-3751324 312902.8, -3779074 328924.4, … WARM_ALB
## 3 WARM 3 ALB ((-3751324 409032.1, -3779074 425053.6, … WARM_ALB
## 4 WARM 5 ALB ((-3723574 264838.2, -3751324 280859.7, … WARM_ALB
## 5 WARM 6 ALB ((-3723574 360967.5, -3751324 376989, -3… WARM_ALB
## 6 NPSW 7 ALB ((-3723574 1418389, -3751324 1434411, -3… NPSW_ALB
Each sf object has the following columns:
province is the longhurst province codecellsID is the unique ID for each hexagonal PUspecies is the species code (e.g. ALB = Albacore)new_feature is the subtyped features code (e.g. WARM-ALB = Albacore in the Western Pacific Warm Pool)geometryTo create climate-smart closures, we prioritized protection of areas with:
We determined these areas using the following metrics:
Relative Climate Exposure (RCE) index
\(\large RCE(yr^{-1}) = \frac{{\text Slope \hspace{2mm}} (^{\circ}Cyr^{-1} \hspace{2mm} 2050 - 2100)}{{\text Seasonal \hspace{2mm} Range} \hspace{2mm} (^{\circ}C \hspace{2mm} 2015 - 2020)}\)
Areas with low RCE index values are areas that will be exposed to low levels of warming relative to their current seasonal fluctuations of temperature. Areas of low RCE index values are prioritized for selection.
Climate velocity
\(\large Climate \hspace{2mm} velocity \hspace{2mm} (kmyr^{-1}) = \frac{{\text Slope \hspace{2mm}} (^{\circ}Cyr^{-1} \hspace{2mm} 2050 - 2100)}{{\text Spatial \hspace{2mm} gradient} \hspace{2mm} (^{\circ}Ckm^{-1} \hspace{2mm} 2050 - 2100)}\)
Areas of slow climate velocity are areas that are more likely to retain their current environmental conditions and consequentially, their biodiversity. Areas of slow climate velocity are prioritized for selection.
The temperatures used to calculate the aforementioned metrics were derived from 12 Global Circulation Models (GCMs) from the Couple Model Intercomparison Project Phase 6 (CMIP6; https://esgf-node.llnl.gov). The GCMs are under these three climate scenarios:
We intersected the RCE values with the PUs using fClimateInt(). The relevant scripts for creating the climate-smart features layer are found in ~/scripts/ with the prefix 07.
# function is found at:
source("scripts/07a_Climate_fClimateInt.R")
# function is ran at:
source("scripts/07b_Climate_fClimateIntRun.R")
We intersected the climate velocity values with the PUs using the same function, fClimateInt().
The cost layer represents the opportunity cost to fishing of establishing the proposed fisheries closures.
The cost layer for each planning unit was represented in US Dollars, and was calculated by multiplying the:
# function is found at:
source("scripts/06a_Cost_fCostPU.R")
# function is ran at:
source("scripts/06b_Cost_fCostPURun.R")
## [1] "minimum: 0.0036497397813946 maximum: 5977.41455078125"
We used log-transformed cost for the figure below:
We join the climate-smart layer with the conservation features layer. Then, we only retain planning units that have RCE and/or climate-velocity values belonging to their lower 25th percentile range. This was how we priortized protection of areas that are predicted to have lower exposure to climate warming and higher retention of biodiversity.
We used the function fFilterQuartile() and the scripts that are relevant for this step are found in ~/scripts/ with the prefix 08.
# function is found at:
source("scripts/08a_Filter_fFilterQuartile.R")
# function is ran at:
source("scripts/08b_Filter_fFilterQuartile25Run.R")
This was done for both the bycatch and large pelagic layers, across different scenarios (SSP1-2.6, SSP2-4.5, SSP5-8.5). Below is an example of how the data looks like, with the following columns:
province: Longhurst provincescellsID: unique ID for each PUspecies: species codenew_features: subtyped features (feature x province)total_area: the total area covered by each subtyped feature (in km\(^{-1}\))velocity: climate velocityvelo_tvalue: transformed values of climate velocity (velocity * 10)RCE: RCE index valuesRCE_tvalue: transformed values of RCE (cube root of RCE)# for bycatch under the climate scenario SSP1-2.6.
bycatchSSP126 <- readRDS("outputs/08_Filter/08b_Filter25/bycatchSSP126_25percentile.rds") %>%
as_tibble() %>%
dplyr::select(-prov_descr, -area_km2, -cellsID.2, -geometry)
head(bycatchSSP126)
## # A tibble: 6 x 9
## province cellsID species new_features total_area velocity velo_tvalue RCE
## <chr> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 WARM 1 Rep-343… WARM_Rep-343… 1664582. 2.16 21.6 0.0740
## 2 WARM 2 Rep-343… WARM_Rep-343… 1664582. 1.37 13.7 0.0626
## 3 WARM 3 Rep-343… WARM_Rep-343… 1664582. 0.531 5.31 0.0613
## 4 WARM 4 Rep-343… WARM_Rep-343… 1664582. 1.33 13.3 0.0765
## 5 WARM 5 Rep-343… WARM_Rep-343… 1664582. 2.91 29.1 0.0689
## 6 WARM 6 Rep-343… WARM_Rep-343… 1664582. 0.743 7.43 0.0585
## # … with 1 more variable: RCE_tvalue <dbl>
# for large pelagics under the cliamte scenario SSP5-8.5.
commercialSSP585 <- readRDS("outputs/08_Filter/08b_Filter25/commercialSSP585_25percentile.rds") %>%
as_tibble() %>%
dplyr::select(-prov_descr, -area_km2, -cellsID.2, -geometry)
head(commercialSSP585)
## # A tibble: 6 x 9
## province cellsID species new_features total_area velocity velo_tvalue RCE
## <chr> <int> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 WARM 1 ALB WARM_ALB 963004. 22.0 220. 3.47
## 2 WARM 3 ALB WARM_ALB 963004. 25.7 257. 2.85
## 3 WARM 5 ALB WARM_ALB 963004. 31.2 312. 3.18
## 4 WARM 14 ALB WARM_ALB 963004. 23.4 234. 3.36
## 5 WARM 36 ALB WARM_ALB 963004. 24.8 248. 3.25
## 6 WARM 56 ALB WARM_ALB 963004. 26.2 262. 3.16
## # … with 1 more variable: RCE_tvalue <dbl>
Aside from these climate-smart closures, we also explored climate-uninformed, wherein the climate-smart features were not included in the analyses. Subsequently, the step of filtering and retention of planning units was also not included.
source("scripts/08c_Filter_fFilterQuartile100Run.R")
This is how the data looks like:
# for bycatch:
bycatch_uninformed <- readRDS("outputs/08_Filter/08c_Filter100/bycatchuninformed_100percentile.rds") %>%
dplyr::select(-prov_descr, -area_km2, -geometry) %>%
as_tibble()
head(bycatch_uninformed)
## # A tibble: 6 x 6
## province cellsID species new_features total_area geometry
## <chr> <int> <chr> <chr> <dbl> <POLYGON [m]>
## 1 WARM 1 ITS-Mam… WARM_ITS-Ma… 912319. ((-3751324 216773.5, -37790…
## 2 WARM 2 ITS-Mam… WARM_ITS-Ma… 912319. ((-3751324 312902.8, -37790…
## 3 WARM 3 ITS-Mam… WARM_ITS-Ma… 912319. ((-3751324 409032.1, -37790…
## 4 WARM 4 ITS-Mam… WARM_ITS-Ma… 912319. ((-3723574 168708.9, -37513…
## 5 WARM 5 ITS-Mam… WARM_ITS-Ma… 912319. ((-3723574 264838.2, -37513…
## 6 WARM 6 ITS-Mam… WARM_ITS-Ma… 912319. ((-3723574 360967.5, -37513…
# for large pelagics:
commercial_uninformed <- readRDS("outputs/08_Filter/08c_Filter100/commercialuninformed_100percentile.rds") %>%
dplyr::select(-prov_descr, -area_km2, -geometry) %>%
as_tibble()
head(commercial_uninformed)
## # A tibble: 6 x 6
## province cellsID species new_features total_area geometry
## <chr> <int> <chr> <chr> <dbl> <POLYGON [m]>
## 1 WARM 1 ALB WARM_ALB 963004. ((-3751324 216773.5, -377907…
## 2 WARM 2 ALB WARM_ALB 963004. ((-3751324 312902.8, -377907…
## 3 WARM 3 ALB WARM_ALB 963004. ((-3751324 409032.1, -377907…
## 4 WARM 5 ALB WARM_ALB 963004. ((-3723574 264838.2, -375132…
## 5 WARM 6 ALB WARM_ALB 963004. ((-3723574 360967.5, -375132…
## 6 NPSW 7 ALB NPSW_ALB 4390870. ((-3723574 1418389, -3751324…
To determine the relative targets per subtyped feature, we checked each one’s IUCN category (IUCN, 2015) using rredlist(). Threatened features (i.e. classified as vulnerable, endangered, threatened) were assigned a fixed maximum target. Features (\(i\)) that were not classified as threatened (e.g. classified as near-threatened, etc.) were assigned a relative target based on their distribution.
Relative target\(_i\) \(= Target_{max} \cdot \frac{PU_i}{PU_{total}} \cdot (Target_{max} - Target_{min})\)
PU_i: planning units of feature \(i\)PU_{total}: total planning unitsTarget_{max}: Maximum target / the fixed maximum target assigned to threatened features (e.g. 1)Target_{min}: Minimum target (e.g. 0.1)Lower targets are assigned to features with broader distribution, and higher targets to features with more restricted distribution.
We reran this step for an array of targets, for both climate-smart and climate-uninformed closures:
| Runs | Maximum | Minimum |
|---|---|---|
| 1 | 1 | 0.1 |
| 2 | 0.9 | 0.1 |
| 3 | 0.8 | 0.1 |
| 4 | 0.7 | 0.1 |
| 5 | 0.6 | 0.1 |
| 6 | 0.5 | 0.1 |
| 7 | 0.4 | 0.1 |
| 8 | 0.3 | 0.1 |
| 9 | 0.2 | 0.1 |
| 10 | 0.1 | 0 |
We used the function fRepresentTarget(), and the scripts are found in ~/scripts/ with the prefix 09.
# the function is found in:
source("scripts/09a_Target_fRepresentTarget.R")
# the runs for climate-smart closures are found in:
source("scripts/09b_Target_fRepresentTargetSmartRun.R")
# the runs for climate-uninformed closures are found in:
source("scripts/09c_Target_fRepresentTargetUninformedRun.R")
For this .rmd, we’ll be showing results using Runs 2 (Max = 0.9, Min = 0.1), 5 (Max = 0.6, Min = 0.1) and 8 (Max = 0.3, Min = 0.1). This workflow shows that it can be replicated using other parameters.
Below is an example of how the data looks like:
target90_commercial <- readRDS("outputs/09_Target/09b-c_TargetRuns/02_Target90/SSP126/target_commercialSSP126.rds") %>%
as_tibble()
head(target90_commercial)
## # A tibble: 6 x 5
## new_features total_area category geometry target
## <chr> <dbl> <chr> <GEOMETRY [m]> <dbl>
## 1 ARCH_ALB 218743. NT POLYGON ((-2086316 -2859363, -2114… 0.185
## 2 ARCH_SKP 18673. LC MULTIPOLYGON (((-2086316 -2859363,… 0.0158
## 3 ARCH_SWO 32011. LC MULTIPOLYGON (((-2086316 -2859363,… 0.0271
## 4 CHIL_YFT 29344. NT MULTIPOLYGON (((8985986 -504196, 8… 0.0248
## 5 non-categ_Long… 424148. NT MULTIPOLYGON (((-476808.4 -2859363… 0.359
## 6 non-categ_Long… 242752. LC MULTIPOLYGON (((-698809.5 -1705812… 0.205
target90_bycatch <- readRDS("outputs/09_Target/09b-c_TargetRuns/02_Target90/SSP126/target_bycatchSSP126.rds") %>%
as_tibble()
head(target90_bycatch)
## # A tibble: 6 x 5
## new_features total_area category geometry target
## <chr> <dbl> <chr> <GEOMETRY [m]> <dbl>
## 1 ARCH_Rep-1347… 88031. VU POLYGON ((-1753314 -2859363, -17810… 90
## 2 ARCH_Rep-2226… 88031. EN POLYGON ((-1753314 -2859363, -17810… 90
## 3 ARCH_Rep-3437… 445489. VU MULTIPOLYGON (((-726559.6 -3003557,… 90
## 4 ARCH_Rep-4127… 96034. CR POLYGON ((-1753314 -2859363, -17810… 90
## 5 ARCH_Rep-5414… 13338. VU POLYGON ((-1864315 -2763234, -18920… 90
## 6 CCAL_Rep-3437… 554861. VU MULTIPOLYGON (((5794721 2139359, 57… 90
new_features: subtyped features (feature x province)total_area: the total area covered by each subtyped feature (in km\(^{-1}\))category: IUCN redlist categorytarget: assigned target in percentageprioritizr (v. 7.0.1.2) uses the following functions to define the problem.
# installing latest version of prioritizr from GitHub
devtools::install_github("prioritizr/prioritizr")
# installing gurobi as the optimizer
install.packages("/Library/gurobi911/mac64/R/gurobi_9.1-1_R_4.0.2.tgz", repos = NULL)
library(priortizr)
library(gurobi)
p1 <- prioritizr::problem(x, feature_list, "cost") %>%
add_min_set_objective() %>%
add_relative_targets(target_list) %>%
add_binary_decisions() %>%
add_gurobi_solver(gap = 0.1, verbose = FALSE)
s1 <- solve(p1)
problem(): we define the datax: data frame having all the data layers (all planning units with their specific ID (cellsID), features data, cost)feature_list: list of the features in xcost: name of the column in x where the cost is foundadd_min_set_objective(): function that allows us to employ the minimum set objective function. This allows prioritizr to create spatial plans that will minimize the cost.add_relative_targets(): target_list list of assigned targets for each featureadd_binary_decisions(): prioritizr will only either select or not select a PUadd_gurobi_solver(): Gurobi was used as the optimizer to solve the problem (at a default 10%)solve(): prioritizr solves the problem p1We used the fPrioritizrRun() function to streamline the runs. fPrioritizrRun() arranges all the data and makes sure that they’re in the right format, then runs problem() and solve(). This function also saves the summaries each plan (cost, % area protected and number of representative features that achieved 30% or more protection).
The relevant scripts for this step is found in ~/scripts/ with the prefix 10.
# the function is in:
source("scrpits/10a_Prioritizr_fPrioritizr.R")
# to run the function with climate-smart data
source("scripts/10b_Prioritizr_fPrioritizrAQMSmartRun.R")
# to run the function with climate-uninformed data
source("scripts/10d_Prioritizr_fPrioritizrAQMUninformedRun.R")
# to create no-regret closures with the climate-smart plans
To create no-regret plans, we used thre function fCreateNoRegret() which intersected the selected PUs in all three climate scenarios. The relevant scripts for this step is found in ~/scripts with the prefix 11.
# the function is in:
source("scripts/11a_NoRegret_fCreateNoRegret.R")
# the runs are found in:
source("scripts/11b_NoRegret_fCreateNoRegretAQMRuns.R")
As mentioned in the previous section, we will only be showing results of Runs 2 (Max = 0.9, Min = 0.1), 5 (Max = 0.6, Min = 0.1) and 8 (Max = 0.3, Min = 0.1).
SSP1-2.6 run
SSP126_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP126.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 90 SSP1.26 267943. 26.0 59
SSP2-4.5 run
SSP245_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP245.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 90 SSP2.45 267309. 26.9 61
SSP5-8.5 run
SSP585_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_SSP585.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 90 SSP5.85 313651. 29.8 60
noregret_target90 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target90.rds")
uninformed_target90 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/02_Target90/solution_uninformed.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 90 uninformed 744951. 73.5 66
SSP1-2.6 run
SSP126_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP126.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 60 SSP1.26 132371. 17.3 54
SSP2-4.5 run
SSP245_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP245.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 60 SSP2.45 131116. 18.0 54
SSP5-8.5 run
SSP585_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_SSP585.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 60 SSP5.85 155441. 19.9 55
noregret_target60 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target60.rds")
uninformed_target60 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/05_Target60/solution_uninformed.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 60 uninformed 369968. 49.0 63
SSP1-2.6 run
SSP126_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP126.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 30 SSP1.26 49676. 8.69 39
SSP2-4.5 run
SSP245_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP245.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 30 SSP2.45 49934. 9.01 39
SSP5-8.5 run
SSP585_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_SSP585.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 30 SSP5.85 57173. 9.96 38
noregret_target60 <- readRDS("outputs/11_NoRegret/11b_AQMRuns/noregretclosures_Target60.rds")
uninformed_target30 <- readRDS("outputs/10_Prioritizr/10b-c_AQMRuns/08_Target30/solution_uninformed.rds")
cost is in USDperc_area_protected is the percentage of the Pacific ABNJ protectedprotected_PU is the number of features that have >= 30% (by area) protection## # A tibble: 1 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 30 uninformed 129125. 24.5 42
Since we ran prioritizr from 10 - 100% maximum target representation (in increments of 10%), we’re showing a summary of how cost and percent protection varied with increasing target representation.
source("scripts/12c_Stat_SomeAnalyses.R")
# log-transformed cost with increasing target representation:
log_plot
# % area protected with increasing target representation:
area_plot
# number of planning units that achieved 30% or more protection (maximum is 66)
protected_features
We also determined Cohen’s Kappa Coefficient between plans created for each run. Below are the summaries for each run as well as a matrix plot of the correlations between plans.
The relevant scripts for this part are found in ~/scripts/ with prefix 12.
# the function for creating the correlation plots is found in:
source("scripts/12d_Stat_fKappaCorrplot.R")
# the runs are found in:
source("scripts/12e_Stat_fKappaCorrplotRun.R")
Summary for Run 2: Maximum Target = 90%
## # A tibble: 4 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 90 SSP1.26 267943. 26.0 59
## 2 90 SSP2.45 267309. 26.9 61
## 3 90 SSP5.85 313651. 29.8 60
## 4 90 uninformed 744951. 73.5 66
Summary for Run 5: Maximum Target = 60%
## # A tibble: 4 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 60 SSP1.26 132371. 17.3 54
## 2 60 SSP2.45 131116. 18.0 54
## 3 60 SSP5.85 155441. 19.9 55
## 4 60 uninformed 369968. 49.0 63
Summary for Run 8: Maximum Target = 30%
## # A tibble: 4 x 5
## target scenario cost perc_area_protected protected_PU
## <dbl> <chr> <dbl> <dbl> <dbl>
## 1 30 SSP1.26 49676. 8.69 39
## 2 30 SSP2.45 49934. 9.01 39
## 3 30 SSP5.85 57173. 9.96 38
## 4 30 uninformed 129125. 24.5 42